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ABSTRACT 

We have examined the soft X-ray plus optical/UV spectrum of the nearby isolated neutron star 
RX J1856 — 3754, comparing with detailed models of a thermally emitting surface. Like previous in- 
vestigators, we find the spectrum is best fit by a two-temperature blackbody model. In addition, our 
simulations constrain the allowed viewing geometry from the observed pulse fraction upper limits. These 
simulations show that RX J1856 — 3754 is very likely to be a normal young pulsar, with the non-thermal 
radio beam missing Earth's line of sight. The SED limits on the model parameter space put a strong 
constraint on the star's M/R. At the measured parallax distance, the allowed range for Mns = 1-5M© 
is Rns = 13.7 ± 0.6 km. Under this interpretation, the EOS is relatively stiff near nuclear density and 
the 'Quark Star' EOS posited in some previous studies is strongly excluded. The data also constrain the 
surface T distribution over the polar cap. 

Subject headings: stars: neutron, equation of state 
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1. INTRODUCTION 

RX J1856 - 3754, discovered by Walter, Wolk, & 
Neuhauser (1996), is the nearest and brightest known neu- 
tron star not showing emission dominated by non-thermal 
magnetospheric processes. As such it offers a unique 
opportunity to study the bare thermal surface emission. 
Measurements of the spectrum can probe the neutron star 
mass (Mns) and radius (Rns), constraining the high den- 
sity equation of state (EOS). Since the discovery, there 
have been several intensive observing campaigns covering 
the optical-UV (HST) and most recently the detailed soft 
X-ray spectrum ( CXO) . An initial 50ks Low-Energy Trans- 
mission Grating (LETG) spectrum showed a broad band 
spectrum remarkably consistent with a simple blackbody 
(Burwitz et al. 2001), although hints of spectral features 
were suggested (van Kerkwijk 2002). A deeper observation 
using 450 ks of Director's Discretionary Time (DDT) was 
made. This unique data set has been the subject of prompt 
study; several authors show that lines in the spectrum are 
undetectable, while pulse searches have placed stringent 
limits on the observed soft X-ray pulse fraction (Ransom, 
Gaensler, & Slane 2002; Drake et al. 2002). These data 
have been variously interpreted, including the widely re- 
ported speculation (based on the X-ray spectrum alone) 
that RX J1856 — 3754 might be a bare quark star (esp. 
Drake et al. 2002). 

Despite the very stringent constraints placed on the X- 
ray pulse fraction, (< 4.5% at 99% confidence, including a 
P search Ransom, Gaensler, & Slane 2002), there is strong 
evidence that RX J 1856 — 3754 is a rotation-powered pul- 
sar, van Kerkwijk & Kulkarni (2001b) discovered an Ha 
nebula surrounding the neutron star, concluding that it 
could be best interpreted as a bow-shock nebula powered 
by a relativistic wind of e^ generated by pulsar spindown. 
The bow shock geometry then provides an estimate of the 
spindown power E = IVtCt = 8 x 10 32 erg/s df 40 (van 
Kerkwijk & Kulkarni 2001b). Adopting magnetic dipole 
braking at constant B, this gives E — 10 34 (i?i2T6) -2 erg/s 



for a surface dipole field 10 12 £?i2 G and characteristic age 
10 6 r 6 y, suggesting B 12 t 6 ~ 3. 

A critical parameter in the discussion of this source is 
the distance, which has been the subject of some contro- 
versy. Initial estimates from HST astrometry gave a paral- 
lax distance of 61pc (Walter 2001). Kaplan, van Kerkwijk, 
& Anderson (2002), however re-analyzed these data, de- 
riving d = 143 pc. A fourth HST observation appears to 
have resolved this discrepancy, giving an overall measure- 
ment of d = 117 ± 12 pc (Walter & Lattimer 2002); we 
adopt this value. 

2. SPECTRAL FITS 

For some time now, it has been clear that the broad 
band spectrum of RX J1856 - 3754 from HST, ROSAT, 
and EUVE data (for a detailed discussion, see Pons et al. 
2002) is inconsistent with a light element (~ Kramer's law 
opacity) H or He atmospheres. H models, for example, 
overpredict the optical/X-ray flux ratio by a factor ~ 100 
(Pavlov et al. 1996). Blackbody, heavy element or compos- 
ite models gave acceptable fits. To produce such ~ Planck- 
ian spectra, one must have nearly isothermal conditions 
at optical depth t u » 1 across the observed band. One 
possibility is that the surface is in a solid or liquid state, 
precluding large temperature gradients through the photo- 
sphere. Theoretical studies to date, while limited, suggest 
that this is not the case unless H is present and B > 10 13 G 
(Lai & Salpeter 1997). In an atmosphere, the easiest way 
to form the spectrum in a small depth range is to invoke 
a rich line spectrum; thus, at first sight, blackbody- like 
spectra should require a high spectral density of opac- 
ity features (lines and edges). This led to the expecta- 
tion that the blackbody-like spectrum of RX J1856 — 3754 
would show many spectral features when examined with 
good S /N LETG spectra. Unfortunately, the initial 50 ks 
exposure already showed that the line features were sub- 
stantially weaker than expected in a simple law-B single 
temperature atmosphere model dominated by heavy ele- 
ments. The new DDT exposure only enhances this conclu- 



2 



Braje & Romani 



sion, placing very strong limits (typical equivalent width 
< 0.02A Drake et al. 2002) on spectral features. 

We will be concerned with the quality of statistical fits 
to various atmosphere models. It is important to note here 
that blackbody spectral fits to the CXO RX J1856 - 3754 
LETG data are, contrary to early reports, not statistically 
perfect fits to a simple Planck spectrum. This conclusion 
was drawn from fits to basic CIAO extractions, which ap- 
preciably underbin the spectrum. Instead we find that 
at a more modest binning (equally spaced ~ 0.7Abins) 
we obtain x 2 /DOF=I.6. As the bin size is increased, the 
X 2 /DOF grows to ~ 4.8, until the number of degrees of 
freedom becomes small. This is a clear signature of spec- 
tral departures on resolved energy scales, and with appro- 
priate binning one indeed finds systematic, grouped resid- 
uals to the Planck function fit at the ~ 10% level. We 
believe that these represent the limit of accuracy in cali- 
bration of the response matrix, as the broad band spectral 
shape is an excellent fit to the Planck function. Drake 
et al. (2002) have reached similar conclusions. Recogniz- 
ing that very subtle departures from a pure blackbody 
may be present in these data, we adopt the conservative 
assumption that these departures are fully accounted for 
by response matrix systematics. 

To accommodate an extended atmosphere, one must 
suppress the spectral features. One possibility is that ex- 
ternal heat sources (such as precipitating magnetosphereic 
e ± ) drive the atmosphere towards isothermality. Sample 
atmospheres showing this effect have been computed in 
Gansicke, Braje, & Romani (2002). A second possibility 
is that the line energies for a given species vary strongly 
across the neutron star surface. For normal pulsar fields, 
B ~ I0 12 G or higher, the strong dependence of the tran- 
sition energies on the local B (eg. Rajagopal, Romani, & 
Miller 1997; Pavlov et al. 1995), coupled with substantial 
> 2 x variation of B across the surface, even for the sim- 
plest dipole models, ensures that such 'magnetic smearing' 
will strongly suppress the phase-averaged line width (Ro- 
mani, in prep.). We discuss briefly here a third possibil- 
ity, that the lines experience variable shifts as the pulsar 
rotates due to Doppler and other dynamic effects. Again, 
the phase-averaged spectrum observed for RXJ1856 — 3754 
would be expected to show broadened and blended lines, 
driving the spectrum towards a Planck curve. 

If, as is required for a normal neutron star, the soft 
X-ray emission of RX J1856 — 3754 is dominated by hot 
polar caps, the rotation of these past the line of sight 
produces phase dependent Doppler shifts (Braje, Romani, 
& Rauch 2000). These are only significant for v sur f ~ 
27tRns/Pns — * c, i.e. Pns < a few ms - Such small Pns 
are not excluded, since the HRC-S wiring error limits ar- 
rival time accuracy and precludes sensitive searches for 
Pns ^ 10 ms. Moreover, the E from the bow shock stand- 
off gives P NS = 4.6(B/10 8 G) 1 / 2 ms, so a low field star, 
having a non-magnetic atmosphere would have a ~ ms 
spin period. 

For a concrete example, we assume Mns = 1.4M , 
Rns = 10 km, and Pns = 1-5 ms (allowed, as the exis- 
tence of PSR 1937 + 21 shows). We have tried both solar 
abundance and iron model atmospheres, tested a range 
of magnetic inclinations a, and computed phase averaged 
spectra using an extension of the Monte Carlo simulation 



code described in Braje, Romani, & Rauch (2000). We 
then tested these models, fitting the most recent CXO 
data, allowing the temperature, observer viewing angle (, 
and interstellar absorption nn to vary. The large L and 
M shell edges in the iron models ensure that these are al- 
ways poor fits. The solar abundance atmospheres have 
a much richer line structure which is more easily blurred 
into a pseudo-continuum. In Figure 1, we display the solar 
abundance millisecond pulsar model, overlaid on the CXO 
data. For comparison the best fit non-rotating model is 
shown in the upper panel. Doppler boosting produces a 
qualitatively acceptable fit below <~ 0.5 keV, but the sim- 
ple blackbody remains statistically an appreciably better 
model (y 2 /DOF=I.6 vs. 3.7 for the binning chosen). We 
must conclude that a simple blackbody fit remains the best 
available, although not yet physically explained. Planck 
emission from a physical neutron star can be compared 
with the data to obtain significant constraints on the neu- 
tron star parameters; we pursue this in the remainder of 
the paper. 
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Fig. 1. — Top: Best fit solar abundance model with no Doppler 
shifts. Bottom: Best fit solar abundance model with Doppler shifts 
for a Pns = 1-5 ms, Rns = 10 km pulsar. 

3. TWO TEMPERATURE BLACKBODY FITS 

A blackbody fit to the CXO data alone results in a tem- 
perature of Too ~ 61 eV with very small statistical un- 
certainty (Drake et al. 2002, and our own analysis). We 
find, as also reported by Drake et al. (2002), that sys- 
tematics (likely in the the effective noted above) 
provide the dominant error. Drake et al. (2002) quote 
61.2 ± 1.0 eV. Taken at face value this with the par- 
allax distance gives a radius as measured at infinity of 
Roo = (1 + z)Rns = 3.8 — 8.2 km. If interpreted as the full 
star radius, this demands exotic equations of state (i.e. 
quark stars). Of course, this is only a lower limit to the 
stellar radius. The optical/UV data points, which closely 
follow a Rayleigh- Jeans spectrum (eg. Pons et al. 2002; 
van Kcrkwijk & Kulkarni 2001a), are most easily inter- 
preted as a second cooler Planck spectrum representing 
flux from the full surface. This has been previously recog- 
nized, but Drake et al. (2002) argue against this interpre- 
tation, citing the absence of the X-ray pulse expected from 
such a hot polar cap/cool surface combination. We have 
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addressed this concern quantitatively, computing detailed 
light curves and spectra. 

3.1. Analytic Two- Temperature Models 

A simple analytic two-temperature blackbody fit delin- 
eates the basic model parameters. For a range of effective 
(cold surface) radii, we fit Ti lot , T co id, the hot area, and 
the absorption column density. In Figure 2, we plot the 
minimum x 2 as a function of cold radius. The optical-UV 
data points fix Rj?j S T co id. The fit becomes poor when T co id 
starts to allow significant Wien peak contribution to the 
CXO X-ray band; this sets the minimum stellar radius. 
Formally, there is a maximum acceptable radius beyond 
which low T co id predicts Wien peak curvature inconsistent 
with the ~Rayleigh-Jeans UV data points. 




Fig. 2. — Solid lino: \ 2 (257 DOF) as a function of cold sphere 
radius. Points: \ 2 values from the polar-cap model fits. 

3.2. More Realistic Two- Temperature Models 

In addition to the phase averaged spectrum, we have a 
limit on the CXO-band pulse fraction. A more realistic 
model is required to address the detailed spectrum and 
pulsations. We adopt a two-temperature model with two 
opposing hot spots (polar caps) at Thot and the remainder 
of the surface at T co id. The caps' orientation (a and () are 
free parameters. We radiate from these surfaces, tracing 
the photons to infinity to form phase resolved spectra and 
light curves. For details of these Monte-Carlo sums see 
Braje, Romani, & Rauch (2000). 

The analytic model results allow some useful simplifica- 
tions. Since Thot is virtually constant over the full accept- 
able Rns range, we fix this value in exploring the rest of 
the parameter space. Further, the X-ray flux amplitude 
allows an initial estimate of the cap half angle A (which 
depends on a and () that ensures that the pulse forma- 
tion is accurate and results in quick convergence to the 
true minimum. 

The fit parameters of interest are Rns, A, T cc jd, C> a , 
and nn- To explore the sensitivity to several parameters, 
we have computed a model grid. We calculate all mod- 
els for a — 5° to a — 90° in five degree steps; ( = 0° 
to C = 90° in five degree steps; and Rns = 12 km to 
Rns = 20km in one kilometer steps. While not strictly a fit 



parameter, we also vary the stellar mass Mns- The quality 
of the spectral fit turns out to be quite insensitive to the 
choice of a and (. In Figure 3, we display a typical spec- 
trum for a model with Mns = 1-4M Q , Rns ^ 14 km. As 
the radius becomes smaller, the optical-UV flux is under- 
predicted and the fit becomes unacceptable. 
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Fig. 3.— Broad band spectral fit to RXJ1856- 3754. Optical/UV 
data points arc drawn from van Kerkwijk & Kulkarni (2001a) and 
Pons et al. (2002). The dotted lines show the unabsorbed hot and 
cold blackbody components. 

For each a and £ we compute \ 2 as a function of Rns- 
The average over angles is shown by the points in Figure 
2. Note that the minimum value is not at xt = 1) a con ~ 
sequence of the aforementioned systematic errors. This 
means that the errors are not Poisson distributed accord- 
ing to the bin counts. To establish confidence levels (CL), 
we must Monte Carlo according to the observed error dis- 
tribution. We have computed the x 2 distribution about 
the best fit model for each (a, () combination, obtaining 
a histogram of x 2 values. These were almost completely 
insensitive to the angles, so we combined all the x 2 dis- 
tributions to obtain confidence levels to the x 2 increases 
associated with variations in Rns- Examinations of the 
differences between the different angles and between inde- 
pendent Monte Carlo runs show that the estimates of 90% 
and 99% confidence level limits on the radius are uncertain 
by no more than ±0.5 km from systematic and computa- 
tional errors in this procedure. 

3.3. Cap Shape Constraints 

One might question whether a simple, circular uniform 
Thot polar cap is merely an adequate approximation to 
the data. We have fit some alternative surface temper- 
ature distributions T(ry), where r\ is the magnetic co- 
latitude, comparing with the best fit circular cap which 
had fcThot = 62.8 eV, A = 21°. The best fit Gaussian 
T(rf) had a peak temperature of fcThot = 74.8eV and width 
u = 19°, but showed an increase of x 2 ~ 80 over the simple 
cap model, sufficient to exclude at the ~ 99% CL. Adding 
a simple linear Thot to T co id ramp to a uniform cap makes 
no discernible difference until the ramp width is twice that 
of the cap. At this point the best fit model has a cap with 
fcThot = 69 eV and A s=a 12°, but the model is excluded at 
the ~ 90% CL. Finally, we fit both a simple Thot oc cos(?7) 
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model and the surface T(rj) distribution of Greenstein & 
Hartke (1983) which is motivated by magnetic anisotropy 
in the thermal conductivity. Both were excluded by the 
CXO data at very high confidence. All fits were at 1.4M Q , 
best fit stellar radii and (a, £) chosen so that the pulse 
fraction is < 5% for the default cap model. 

Evidently, the CXO data require a quite uniform distri- 
bution of the high temperature excess, and suggests that 
it is induced by exterior heating rather than interior con- 
ductivity. This cap size is substantially larger than the 
A ~ 3°(100ms/PNs) 1 ^ 2 expected for a dipole surface cap, 
unless the period is very small. Higher magnetic multi- 
poles would generally have even smaller open zone caps. 
The large, uniformly heated area is puzzling in the context 
of pulsar surface acceleration models, but might be most 
easily accommodated in the more modern pictures of a 
GR-induced potential that is relatively uniform across the 
polar cap and forms a pair formation front at relatively 
high altitudes (Harding & Muslimov 1998). Perhaps a 
more plausible interpretation invokes a high-altitude ac- 
celeration zone, with the inward-directed 7-rays pair con- 
verting in the closed zone above the polar cap (Wang et al. 
1998), which should give a ~ large zone of uniform surface 
heating. 

4. PULSE FRACTION CONSTRAINTS 

We have seen that the spectral fits are quite insensitive 
to the cap viewing angles. They place a firm minimum on 
the allowed Rns but allow radii that are implausibly large. 
However, the observed pulse fraction depends strongly on 
viewing geometry and, through gravitational focusing, the 
value of Mns/R-ns- In general, we expect a strong X-ray 
pulse when the magnetic axis passes close to Earth's line- 
of-sight (small \a — (\). Of course, for an aligned rotator 
a s=a the pulse can be very small and for a ~ tt/2 there 
is an appreciable region where the pulse is weak as viewed 
from Earth. We have examined the CXO-band pulse pro- 
file of our model grid to find the allowed region of (a, £) pa- 
rameter space in which the model pulse fraction is weaker 
than that observed (< 4.5%); this is shown in the inset of 
Figure 4. 

This exercise is repeated for each Mns, Rns! the allowed 
phase space is larger for small Rns /Mns as gravitational 
bending dilutes the observed pulse. By demanding a cer- 
tain minimum probability that a pulse fraction as low as 
that observed is seen at Earth, we obtain a maximum ac- 
ceptable radius for the neutron star. One important addi- 
tional constraint can be invoked. The non-detection of this 
neutron star in the radio band (Brazier & Johnston 1999) 
puts a very strong bound on pulsar emission directed to- 
wards Earth. Assuming that this is a normal radio pulsar 
(consistent with E inferred from the bow shock), we must 
conclude o priori (independent of the X-ray data) that our 
line of sight lies outside of the pulsar radio beam and that 
\a — C| is not small. The region excluded depends on the 
size of the radio beam (O ra d), which in turn depends on 
the spin period. A typical estimate is (eg. Rankin 1993) 

O rad - 5.8°(Pns/1s)- 1/2 (1) 

at the radio frequency v = 1 GHz. With Pns = 0.3 s, this 
gives O rad ~ 10.6°. Radio limits on RX J1856 - 3754 arc 
strong at even lower v where radius-to- frequency mapping 



produces larger O ra d, and for this nearby object, the radio 
luminosity constraints are so severe that we are unlikely 
to intersect even the faint fringe of the radio beam. Both 
effects argue for O ra( j larger than that above. Accordingly, 
in the Figure 4 inset we show by diagonal lines the regions 
near a = ( excluded by the lack of radio detection for 1 x 
and 2x this fiducial beam size. Shorter periods allow an 
a priori exclusion of even more phase space. The point 
is that from the lack of radio detection we should have 
expected a low X-ray pulse fraction. The probability for 
obtaining pulses as weak as observed is then set by the 
fraction of the remaining allowed solid angle. These frac- 
tions arc plotted in Figure 4 with and without the radio 
prior. The allowed radius range depends on mass and in 
Table 1, we give the 90% and 99% CL upper bound on the 
neutron star radius for these priors and several neutron 
star masses. 
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Fig. 4. — Fraction of the sky allowed by pulse fraction constraints 
as a function of stellar radius. The solid line assumes no priors; 
the dotted line for a prior 1 X rat j = 10.6°; and the dashed line 
for 2 X © ra( j. The inset shows the allowed (a,C) parameter space 
shaded in light, medium, and dark gray for 12, 14, 16km, respectively 
at Mns = 1-4Mq. The lines in the inset depict the parameter space 
excluded by radio prior. 

5. EQUATION OF STATE CONSTRAINT 

We see that the X-ray/optical data, using the spectral 
and pulse fraction arguments, give a range of allowed stel- 
lar radii for each mass. Strictly speaking, the minimum 
and maximum radius CL have somewhat different inter- 
pretations, but it is interesting to place these bounds in 
the Mns — Rns plane to compare with the predictions of 
various EOS. 

In Figure 5, we show the combined constraints, assum- 
ing O ra d = 21.2°. The spectral radius lower bounds at 90% 
and 99% CL approximate curves of constant Mns/R-ns- 
The pulse fraction upper bounds (90%, 95%) rapidly drive 
one to small radii at low masses. For Mns S L3M q no si- 
multaneous solutions are allowed consistent with the 90% 
bounds. At Mns = L5M© the allowed range from the 
fit is quite small, Rns = 13.7 ± 0.6 km. The additional 
uncertainty in the distance actually dominates the errors 
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(arrowed bar). 

For comparison, several EOS curves (after Lattimer & 
Prakash 2001) are plotted. We see that large radius (stiff 
at nuclear density) EOS are preferred. Formally, the rel- 
ativistic field theoretical model by Miiller & Serot (1996) 
and the model GS2 by Glendenning & Schaffner-Biclich 
(1999) arc the only modern models allowed (the origi- 
nal PS model of Pandharipande & Smith (1975) is also 
allowed). We note that the GS models are very sen- 
sitive to the K meson potential; for example, GS1 is 
strongly excluded. Interestingly, no potential or varia- 
tional method computations agree with the formal over- 
lap. If one includes the distance uncertainty, a few more 
intermediate radius models are not excluded at the 90% 
CL. However, even including the distance uncertainties, 
all quark star models are excluded at the ~ 95% level for 
Mns < 1.5M Q and are only barely consistent at the high- 
est allowed masses. Improved HST parallax measurements 
could boost this exclusion to the 3cr level. 
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Fig. 5. — Radius constraints for different possible RXJ1856 — 3754 
masses. The triangular gray shaded region represents the formal 
90% CL overlap from the pulse fraction and spectral constraints. 
The two-sided arrow represents the systematic range induced by 
the distance uncertainties. Equation of state curves and labels are 
drawn from Lattimer & Prakash (2001). See this paper for EOS 
labels and references. 



The model GS2 differs from the MS models in that it in- 
cludes a kaon condensate in the core. One principal effect 
of exotic interior condensates is to enhance neutrino cool- 
ing after ~ 100 y. This can also be achieved without exot- 
ica for a proton fraction > 11% through the direct URCA 
process (Lattimer et al. 1991). Our inferred T co id, inter- 
preted as the signature of the cooling from the initial heat 
of formation, corresponds tolw 3.2 x 10 31 erg/s for our 
best fit radius. Such luminosities are reached in ~ 5 x 10 5 y 
(the preferred age of Walter & Lattimer 2002) when en- 
hanced neutrino cooling can occur, but are achieved after 
~ 1.3 x 10 6 y for stars with low density (stiff) cores Tsu- 
ruta et al. (2002). Thus, depending on the actual stellar 
age, this cool surface may be seen as weak evidence for 
an exotic composition and/or significant softening of the 
EOS at very high densities. 

The prospects for further tests of the ideas in this paper 
hinge on the detection of a pulse from RX J1856 — 3754 
along with measurement of the pulse fraction and the pe- 
riod derivative. Given that significant allowed pulse frac- 
tion parameter space lies just slightly below the CXO 
LETG detection threshold, the prospects for a pulse mea- 
surement with the recently completed 58 ks XMM obser- 
vation are quite good. Even if only thermal, the phase 
resolved spectrum should provide important constraints 
on the cap temperature, size and orientation. In partic- 
ular our thermal model predicts that the pulse fraction 
should increase by 70% from 0.15 keV to 0.5 keV, where 
the XMM data should still deliver ~ 0.8 PN camera cps, 
allowing detection of pulse fractions well under 1%. 

We are grateful to Boris Gansicke for collaboration on 
the atmosphere models used in §2; and to Herman Mar- 
shall for sharing an independent LETG response matrix. 
This work was supported in part by NASA grant SP2- 
2002X. 
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Table 1 

Pulse Fraction 90%/99% Confidence Limits 



M NS 


No Prior 


e rad = 10.6° 


©rad — 21.2° 


(M ) 


(km) 


(km) 


(km) 


1.2 


10.3/14.6 


11.3/15.6 


12.0/16.3 


1.3 


10.5/16.7 


11.7/17.6 


12.6/19.2 


1.4 


12.2/17.8 


13.1/18.9 


13.6/19.9 


1.5 


13.4/19.4 


13.8/21.0 


14.3/21.3 


1.6 


14.0/20.0 


14.6/21.4 


14.9/21.2 



